NASA Contractor Report 178207 
ICASE Report No. 86-76 


ICASE 


A NUMERICAL ALGORITHM FOR OPTIMAL FEEDBACK GAINS 
IN HIGH DIMENSIONAL LQR PROBLEMS 


H , T . Banks 
K. Ito 

(NASA-CE-178207) A NUMERICAL ALGORITHM FOR N87-14C60 

OFIIMAL FEEDBACK GAINS IN HIGH DIMENSIONAL 
ICfi PECBIEMS Fiiial Report (NASA) 32 p 

CSCL 72B Onclas 

G3/64 43956 

Contract Nos. NASl-17070, NASl—18107 
November 1986 

institute for Computer Applications in Science and Engineering 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


IVIASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 


A Numerical Algorithm for Optimal Feedback Gains 
in High Dimensional LQR Problems 


H.T. Banks 
K.Ito 


ABSTRACT 

We propose a hybrid method for computing the feedback gains in linear 
quadratic regulator problems. The method, which combines use of a Chandrasekhar 
type system with an iteration of the Newton-Kleinman form with variable 
acceleration parameter Smith schemes, is formulated so as to efficiently compute 
directly the feedback gains rather than solutions of an associated Riccati equation. 
The hybrid method is particularly appropriate when used with large dimensional 
systems such as those arising in approximating infinite dimensional (distributed 
parameter) control systems (e.g., those governed by delay-differential and partial 
differential equations). Computational advantages of our proposed algorithm over the 
standard eigenvector (Potter, Laub-Schur) based techniques are discussed and 
numerical evidence of the efficacy of our ideas presented. 


Research supported in part by the National Aeronautics and Space Administration 
under NASA grant NAG-1-517, the Air Force Office of Scientific Research under 
contracts No. AFOSR-84-0398 and AFOSR-85-0303, and the National Science 
Foundation under NSF grant MCS-8504316. Parts of the research were carried out 
while the authors were visiting scientists at the Institute for Computer Applications 
in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 
23665, which is operated under NASA contracts NASl-17070 and NASl-18107. 


i 



- 1 - 


1. Introduction 


A great deal of effort in recent years in control of distributed systems has 

focused on approximation techniques (for example, see [l]-[6], [8], [11]-[14], [16], [18], 

[21], [22], [24], [26], [32]) to reduce inherently infinite dimensional problems to (large) 
finite dimensional analogues. Relatively little attention has been given to the 

development of efficient computational methods for the resulting large but finite 
dimensional control problems. In this paper we consider such questions for one 
classical formulation of the feedback control problem, the well-known linear quadratic 
regulator (LQR) problem. 

There are several approaches one can take in such an endeavor. With the 
emergence of new computer architectures (vector and parallel), one exciting possibility 
involves the development of new algorithms to be used with nonsequential computers. 
While we are currently investigating ideas in this direction, our presentation here is 
to report on some of our efforts to develop better algorithms for use with 

conventional serial computers. 

As is well-known, the LQR problem can be reduced to the solution of a 
matrix Riccati equation in order to construct the feedback gain matrix. The most 
widely available method for solution of the Riccati equation is the Potter method 
[30], which is based on obtaining the eigenvectors and eigenvalues of an associated 
2nx2n Hamiltonian system when the underlying dynamical control system is of 
dimension n. A related, but improved version involving computation of Schur vectors 
for the system was proposed by Laub in [27]. While both of these "eigenvector" 
methods can be used satisfactorily (for a discussion of real advantages offered by the 
Laub-Schur approach over Potter’s method, see [27]) for systems with n relatively 
small, say n < 100, the computational effort (and time) grows like n* and becomes 
prohibitive for large systems. More recently, the idea of using Chandrasekhar systems 
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([9], [20], [33, p.304-310], [36]) when the number of states is large compared to the 
number of control inputs (exactly the situation in a number of cases where one 
approximates a distributed system) has been suggested by a number of authors [33, 
p.309], [7], [8], [17], [31]. However, as we shall discuss below, there can be numerical 
difficulties in using the Chandrasekhar approach. On the other hand, it is known 

that iterative methods such as the Newton algorithm as formulated by Kleinman in 
[23] can be quite efficiently implemented (even for some large systems) if good 
initial estimates are provided and if one can solve efficiently the resulting Lyapunov 
equations. In this presentation, we discuss the formulation and numerical testing of a 
hybrid method that represents an attempt to combine the good features of the 
Chandrasekhar approach (growth like n in computational effort) with those of the 

Newton-Kleinman (quadratic convergence when good initial estimates are provided) 
along with innovative use of the Smith algorithm for solution of Lyapunov 
equations. 

We expect these ideas to be quite useful in design of control laws for some 
of the models currently being investigated in connection with large flexible structures 
as well as in some of the population dispersal and control studies that we are 
currently pursuing with biologists and ecologists. Some of the large flexible structures 
involve rather sophisticated distributed parameter models (e.g. see [4], [34]), especially 
when one wishes to include complicated damping mechanisms involving time or 

spatially related hysteresis [6], [34] or nonlinear effects [19]. For such models, the 
computational tasks can be rather demanding whether one is carrying out parameter 
identification [3] or feedback control calculations with traditional eigenvector based 
methods (the authors of [6] have indicated experiences with runs requiring 9 hours of 
VAX time when using an approximate system with dimension equivalent to n = 

238 ). 

For our presentation, we assume that one has used their favorite 
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approximation scheme (finite-elements, spline, spectral, etc.) to reduce the problems of 
interest to an LQR problem with finite dimensional system. More precisely, 
throughout our discussions we consider the LQR problem; minimize the cost 
functional 


(l.I) 


J(u) 


“|Cx(t)|2 + |u(t)f)dt 
0 


subject to the state dynamics 

x(t) = Ax(t) -I- Bu(t) , x(0) = Xq . 

Here A e R"’^", B e and C e RP^". (We have, without loss of generality, for 

our discussions here normalized our problem so that the control term in the cost 

functional (1.1) appears with a weighting matrix I.) We shall assume that (A,B) is 

stabilizable and (C,A) is detectable [25], [37]. Then the optimal feedback control for 
the LQR problem involving (1.1) is given by 
u(t) = -B'^Px(t) 

where P is the unique non-negative symmetric solution of the algebraic Riccati 
equation 

(1.2) A"^P + PA - PBB'^P + C'^C = 0. 

In this paper we propose an algorithm which leads to direct calculation of the 
feedback gain matrix K = B"^P without computation of P. In addition to providing 

substantial savings in computational time over eigenvector methods, our algorithm 
requires much less storage and can easily be implemented to take advantage of 
special structures (e.g. sparsity) in the system matrix A. 

To outline the steps in our algorithms, we first recall that the optimal 
feedback gain K is given by the limit K = lim K(t) as t of solutions of the 

Chandrasekhar system [9], [20] 
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K(t) = -BTL^(t)L(t) , K(0) = 0 


(1.3) 

L(t) = -L(t)(A-BK(t)) , L(0) = C, 

where K € and L € RP""". In fact (see [20], [37]) P = lim J° L'^(s)L(s)ds as t ^ 

The first step in the proposed hybrid algorithm involves a numerical integration 
of (1.3) backward in time on an appropriate interval [-t^O]. For the second step, the 
value K(-tf) obtained through this numerical integration of the Chandrasekhar system 
is then refined by use of the Newton-Kleinman algorithm [23], if we use K(-tj.) as an 
initial value Kq for the Newton method. 

To motivate our effort in these two steps, we note that the convergence 
K(t) -♦ K as t can be very slow when the eigenvalues of A-BK lie close to the 

imaginary axis. Moreover, L=0, Koo arbitrary are solutions in the asymptotic limit 

sense to (1.3). That is, if we denote by f(K,L) the right side of system (1.3), then K,*, 
arbitrary and L=0 are solutions of f(Keo,L) = 0. Hence K(t) -► Ka,, L(t) -* 0 as t 
doesn’t, in general, have a unique limit numerically . Thus, as is pointed out in [33,p. 
316-318], if one is to use the Chandrasekhar approach alone, one needs a very 
accurate numerical solver for (1.3). This can be computationally quite expensive if we 
are dealing with a large system and/or a stiff system. Hence, we propose to use a 
rather crude, fast integration method for the Chandrasekhar component of our 
algorithm and take the resulting numerical solution K(-tf) as a start-up value for the 
Newton iterations. If this crude estimate from the Chandrasekhar step is a 
sufficiently good initial guess, then we can expect to meet the Newton-Kleinman 
requirements that A-BKq be a stability matrix and to obtain quadratic convergence in 
this second component of the algorithm. 

The first step of our hybrid method requires the solution of n(m+p) 
simultaneous equations, while each iteration of the usual Newton-Kleinman step 
requires the solution of a Lyapunov equation for the nxn symmetric estimates of P. 
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However, as we shall see below, one can use factorization ideas [20] and the Smith 
method [35] for Lyapunov equations to reformulate the Newton-Kleinman method as 
a direct iterative method for the mxn gain K, thereby providing additional 
computational advantages. To speed up our calculations and improve convergence in 
the Smith algorithm, we propose a variable stepsize Smith method to solve the 

Lyapunov equations as described in Section 4 below. In Section 2 we outline a 

numerical scheme for the Chandrasekhar system, while the reformulated 

Newton-Kleinman iterative procedure to compute directly the gain K is detailed in 

Section 3. Finally, in Section 5, we discuss further some advantages and disadvantages 
of the proposed algorithm and report on our experience with several numerical 
examples to illustrate the feasibility of our hybrid approach. 
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2. Numerical Solution of the Chandrasekhar System 

Wc return to consider more closely the Chandrasekhar system: 

(2.1) k(t) = -B'^L'^(t)L(t) , K(0) = 0 

(2.2) L(t) = -L(t)(A - BK(t)) , L(0) = C , 

where A e L,C e and K,B^ 6 We first observe that the second 

equation (2.2) is linear in L. In cases where A arises from a discretization or 
approximation of partial differential equations, equation (2.2) tends to be a stiff 
system and thus it is advisable to use an implicit numerical scheme. We propose here 
the second order Adams-Moulton algorithm [10, p.235]. A second observation is that 
the right hand side of equation (2.1) is independent of K and thus an explicit 
scheme is appropriate; we propose the second order Adams-Bashforth algorithm [10, 

p.226]. 

These observations lead us to propose the following algorithm for the 
Chandrasekhar system (2.1), (2.2): Given a step size h>0, approximations Kj and Lj to 
K(-ih) and L(-ih) are generated by 

(2.3) = K; + h(|BTL^ L. - ^bTl^^Lj.,) 

( 0 ) ( 0 ) 

(2.4) = (Kf^; + K.)/2 

(2.5) - Lj + + L,)(A.BKj°’^ 

(2.6) = K. + 

where Kq = 0 and = Lq = C. 

Several remarks may be useful at this point. 


( Remark H The stiffness of the matrix A dictates the choice of stepsize h. 
( Remark 2) The predicted values and the corrected values satisfy 


-7- 


- K 


( 0 ) 


"i+1 


i+1 


= ll bT _d 
2 dt 


2 

(L^L) 


Ti 


- 2L?L, + Lj,L,.i) 


and this relationship can be used for stepsize refinement, i.e., to give local bounds 
depending on stepsize which can be used in error control. 

( Remark 31 The formula (2.5) can be rewritten as 


(2-7) L.„ . L,(I + h A.)(i . 

. 2L,(I - - Lj 

where Aj = A - . Defining H = I - ^A we have that I - ^A| = H + 

where B e Rnxm ^ . Thus by the Sherman-Morrison-Woodbury 

formula [29, p.50] (used frequently when updating an nxn matrix by rank m 
matrices) 


(2.8) (I - = H-l - ^-iB(I + 

where I + ^k|°|^-1b € R'"’'™. The matrix in (2.8) can be computed by 


(2.9) 


- KiH-> + ^IbT'lT'LP'* - ^TlT,L|.,H-') 
K,H-> t + BVLiH-‘). 


Hence we see from (2.7) - (2.9) that the step (2.5) only involves the operation LjH"^ 
plus inversion of an mxm matrix I + H'^B. Thus the step (2.5) can be 

reformulated so that it requires only an mxm matrix inversion plus matrix-vector 

t- 

multiplications if the LU decomposition of H = I - £A is computed a priori. This 
procedure can be most advantageous computationally when m and p are small 
compared to n. 

( Remark 4) For some problems one might wish to use a completely implicit scheme 
in place of (2.3) - (2.6) to enhance stability and reduce sensitivity to step size choice. 
Then one might consider iterations generated by 
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L: 


HI = Li + + L;)(A-BK;j:i^ 


( 2 . 10 ) 


^ i+1 ^+1 ^ i+1 ^ ® 


*^“(5 = ( K ,« + K ,)/2 


and thus produce iterates with limits (as j ~* ") Kj^^, 


( 2 . 11 ) 


Li+i - L, + ^Li+, + L,)(A.B(K,^., + 

K.,„ - K, + 


satisfying 

Kj)/2) 
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3. An Iterative Method for Computing the Optimal Feedback Gain K 

A widely used iterative method for finding the non-negative solution of the 
algebraic Riccati equation (1.2) is Newton’s method as modified by Kleinman [23]. We 
show that this method can be reformulated so that, when combined with a factored 
form of the well known Smith method [35], one can compute directly a sequence of 
iterates Kj for the feedback gain K. 

First we recall the Newton iterative algorithm as formulated by Kleinman: 

(1) Choose a gain matrix Kq so that A-BKq is a stability matrix; 

(2) Update Kj by = B^Pj where Pj is the solution of the Lyapunov equation 

(A-BK;)'^Pj + Pj(A-BKj) + kTKj + C'^C = 0. 

It can be shown that 0 5 Pj^^ « Pj for any i, and P = lim Pj where the 

convergence is quadratic. This algorithm can be viewed as an iterative method for 

the gain K, i.e. K = lim Kj where = F(Kj) with F(K) = B'^P and P is the 

solution of the Lyapunov equation: 

(3.1) (A-BK)"^P + P(A-BK) + K% + C'^C = 0. 

Thus, in order to calculate F(K) one has to solve the Lyapunov equation (3.1) for the 
symmetric matrix P. However, one can form an alternative version that allows one to 
directly calculate F(K) using the Smith method for a Lyapunov equation in X of the 
form 

(3.2) S'^X + XS + D'^D = 0 

where S € is a stability matrix and D € 

To this end, we replace step (2) in the Newton-Kleinman method by: 

(2') For i ^ 1, update Kj by = Kj - B'^Zj where Zj = Pj_^ - Pj is the solution of 

the Lyapunov equation 

(3.3) (A-BKj)^Zi + Zj(A-BK;) + dTDj = 0 
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with D. = K. - Kj.^. 

The method with (2') offers several advantages over that using (2). The 
Lyapunov equation in (2') has fewer inhomogeneous terms than does the one in (2) 
and the term Dj has rank m which depends only on the number of inputs (controls) 
to the system. In the proposed modified Smith method described below, one is able 
to compute directly the mxn update matrix J‘ = B'^Z^ without computing Zj (see (3.12) 
below). Since Zj 0 as i -* <» is expected, choosing the start-up value Jq = 0 in the 
factored Smith algorithm (where J' = B"^Zj is computed as the limit as k " of a 

sequence JJ^) is a natural as well as convenient choice. 

Note that the step (2') requires that one have Z^ = Pq - in hand and 

hence we must start this procedure with Pq, Pj^ (and Kq,K^) given whereas (2) requires 

only that one start with Kq given. Then is computed by Kj^ = B'^Pq with Pq the 
solution of 

(A-BKq)TPq + Pq(A- BKq) + kJ Kq + C'^C = 0. 

Since our Smith algorithm below is formulated to solve Lyapunov equations of the 
form (3.2), we can, to maintain this form, initially solve the equation twice. That is, 
if we solve for Zq the solution of 

(3.4) (A-BKq)T^ + Zo(A-BKq) + kJ Kq = 0 
and for Zq the solution of 

^ A A fp 

(3.5) (A-BKq) Zq + Zq(A-BKq) + C C = 0 , 

then we can obtain by = B"^Zq + B^^Zq . Since the Smith method as 
formulated here actually returns B^X where X is the solution to (3.2), we thus will 
use this Smith algorithm twice (with S = A-BKq), once with D”^D = KqKq, once with 
0*^0 = C’^C and then simply add the solutions to obtain K^. 

We turn next to the desired factored form of the Smith method as applied 
to equation (3.2). Let Xq be an arbitrary nxn symmetric matrix and let a sequence 
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{Xj^} of nxn symmetric matrices be generated by 

(3.6) = UjXi^U, * Y, , 

where r is a postive constant (the Smith stepsize) and 


(3.7) Ur = (I-rS)'' (I+rS) 

(3.8) Yr = 2r(I-rSV^ D'TD(I-rS)'^ . 

Then one can argue that {X^} converges to X, the solution of (3.2). The method and 
its analysis is based on the observation that for any positive constant r, the equation 
(3.2) is equivalent to 

X = Uj XUr + Yr 

which can be used to define a contraction map in the obvious manner. 

We modify this standard formulation of the Smith method to suit our 
particular needs here (computing B^X instead of X). From (3.6) we have 
X|,« ■ = u7 (X^ - X^.,)U, . k s I. 

Hence, if Xj^ - X^_j = have a factorable difference), then 

Xk+, - X^ = U7 Mf M^U, - (M^U/(M,U,). 

If the start-up value Xq is zero, then we can write 

(3.9) Xj - Xq = 2rM^Mj , = D(I-rS)‘^ . 

By induction on k, (3.6) is then equivalent to 


(3.10) . M|.U, 

(3.11) X^,., = X^ t . 

In this manner B^X can be obtained as the limit of = B”^X^ where is generated 
by 

Jk+i = Jk + 2rBTMT^iM^+i . 

Thus, the update step (2') is carried out by the following Smith algorithm: 


(3.12) (/) 

(3.12) (//) 


Set Sj = A-BKj and D = Kj - 

Choose a positive constant r and form Uj. and by (3.7) and (3.9) 
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(3.12)(m) 


with S = Sj ; put Jq = 0 and Ji = Jq + 
Iterate for k = 1, 2, .... on 


Mk+i = 

^k+l “ ^*'®^^k+l^k+l • 

In summary, we have described in this section a Newton-Kleinman scheme 
combined with the Smith method for the resulting Lyapunov equation at each step in 
the Newton-Kleinman. We have reformulated the Newton-Kleinman iteration and 
factored the Smith algorithm so as to result in algebraic savings in computing 
directly the gain estimates Kj. 
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3. The Smith Method and Variable Steosize 


As is well-known, the rate of convergence in the Smith method discussed in 
the last section depends upon the choice of the acceleration or step parameter (Sec 
[33, p.291-297] for several discussions. Note that our parameter r is the negative 
reciprocal of the parameter in Russell’s discussions). To increase speed in 

convergence, one may employ the accelerated Smith method [33], [35] which can yield 
quadratic convergence as compared to the linear convergence obtained with (3.6). 
However, unlike (3.6), the accelerated Smith method is not self-correcting [33] and 

here we propose to speed up convergence in an alternative way which has proved 

both reliable and efficient in some of our numerical tests. Specifically, we propose to 
use a succession of acceleration parameter values rj (much in the spirit of other 
well-known iterative methods such as alternating directions [15], [28]) to accelerate 

convergence in the basic Smith method. Our formulation of this "variable stepsize" 
Smith method is based upon the observation that for fixed r > 0 and k > 1, the 
Smith algorithm can be written as 

S“^X^ + X. S + D“^D = E. 

(4.1) 

s (I + rS)'^Mj M^(I + rS) 

where M,^ is defined by (3.9) and (3.10). To see this, we note that from (3.6) we 
have 

(I - rS)'^X^(I - rS) = (I + rS)’'X^.j(I + rS) + 2rD'^D , 
or 

(I + rS'^)X^(I + rS) - (I - rS)'^X^(I - rS) + 2rDTD 
= (I + rS)T(Xfc - X^.i)(I + rS) . 

Hence, from (3.11) we obtain 

2r(S'*'X^ + X^S + D'^D) = 2r(I + rS)'^MjMk(I + rS) 
which implies (4.1). Moreover, from (3.9) and (3.10) we have 
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(4.2) Ej^ = (Uj )kDTD(U^)>' , k ^ 1 . 

Thus, if we use the iteration (3.6) with acceleration parameter r^ for iterates, we 
obtain an iterate and equation error satisfying 

+ D'^D = E^^) 

(4.3) . k 

eW = (UJ) 1 D'^D(U, ) ^ 

Let us define the difference = X - X^^^ where X is the sought-after solution of 
(3.2). Then it is readily seen that satisfies a Lyapunov equation similar to that 
of (3.2) : 

(4.4) S^E + ES + e(^) = 0. 

If we next apply the iteration (3.6) kj times with acceleration parameter to the 
residual equation (4.4) we obtain 

(4.5) S'^X^2) + + e(^) = e(2) 

where X^^^ is the final iterate using and the equation error E^^^ is given by 
e(2) = (U^)‘'2 e(^)(U )‘'2 . 

If we proceed to define the difference E^^^ = X - (X^^^ + X^^^), then from (4.4) and 

(4.5) we see that E^^^ satisfies a Lyapunov equation 

S'^E + ES + e(2) = 0 . 

We continue this procedure, using a sequence of acceleration values (rj along 
with corresponding iteration counts {kj to produce a sequence {X^‘^} of nonnegative, 
symmetric matrices. For i > 1, we have 

(4.6) S'^X(‘) + X^S + = e(‘) 

(4.7) e(‘) = (U^)’'‘ e(‘-^)(Uj,)''* , e(°) = D'^D . 

Thus, if Xj s i; X^‘^, then X. satisfies 

i=i 

(4.8) S'^Xj + XjS + D'^D = E^j) 
and hence X. i X and Xj_^ « Xj , j ? 1. 

Using arguments similar to those in [33, p.291-297] one can show that for 
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0 < r 5 r ^ R, with r, R positive constants, there exists a constant u, 0 < u < 1, 
depending only on r and R, such that for u < p < 1, 

|Uj? I « M(p)p>' , k = 0, 1, . . . , 

where M(p) is independent of r. Thus, if r 5 rj ^ R, then for any 0 < 6 < 1, there 
exists an integer k(6) such that for kj ^ k(S) , i ^ 1, 

(4.9) |uj‘| « 1 - 6. 

Hence, using (4.7) we have 

|E(j)| « (l-£)^j |D|2 

so that 0 as j ~ and therefore Xj -» X as j 

For the hybrid method proposed in this paper, we have combined the 
variable stepsize method just outlined with the reformulated Smith method of (3.12). 
We then obtain the following algorithm for solving for the feedback gains K^; 


Algorithm M.lOk 

Set Sj = A - BKj, = Kj - and Jq = 0. For given acceleration parameters r^, 

rj, . . . , and iteration count indices k^, kj, . . . , we iterate on j = 1, 2, . . . , in 
the following steps: 

(4.10a) Compute Uj. = (I - rjSj)‘^(I + rjS;) and 

M, = D.(I - , 

= J 0 + 2rjBTM[Mi 

(4.10b) Iterate for k = 1, 2, . . . , k j - 1 in 

Mk+l - M^U,. 

^k+1 ~ ^k+l^k+1 

Compute (I + r^Sj), set Jq = Jj^ and return to (4.10a) with 

j = j + 1 . 


(4.10c) 
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The steps in Algorithm (4.10) are to be repeated, i.e. iteration through 
rj^, rj, . . . , until some convergence criterion is met. In performing the steps in 
(4.10b), one can use the Sherman-Morrison-Woodbury formula in a procedure such as 
that outlined in (2.7) - (2.9) in section 2. 

As is often the case in variable stepsize algorithms, the choice of the 
acceleration parameters rj^, rj, . . . , and the associated iteration counts k^, kj, . . . , 
provides both freedom and some frustration in the search for "best" choices. If one 
follows the guide provided by ADI methods (see [15, p.37]), one might choose a set 
of values rj to be used in some cyclic order. The best choices of values for the rj 
often depend on the eigenvalues of Sj = A - BKj. For example, consider the case 
where Sj has only real eigenvalues Xj, each with multiplicity mj, j = 1, 2, . . . , m. 
Then a choice of tj = -1/Xj and kj = mj in the algorithm produces convergence in a 
finite number (m) of steps. That is, this choice yields = 0 in (4.8). 

Of course, the complete eigenstructure of Sj will not be known (nor do we 
suggest that any sophisticated analysis along these lines be included with each use of 
Algorithm (4.10) to obtain the gains Kj). A possible alternative is to use one of the 
polynomial acceleration methods [15, Chp 3, 4]. 

In closing this section, we note that the analogies of our variable step Smith 
method with the ADI methods used to solve partial differential equations can be 
made a little more precise. Briefly, in ADI spliting methods [28, p.146-148], one 
attempts to solve a discretization of the evolution equation 

$ = A4> + f 

when A > 0 can be written A = A^ + Aj with Aj ^ 0 (for example, factored into 
components corresponding to spatial discretizations in the x and y directions 
respectively for an equation in a two dimensional spatial domain). This can be shown 
[28, p.l50] to lead to an iterative scheme 



-17- 


(4.11) (I + jAj)(I + ^A2)4>j+^ = (I - jAj)(I - + hfj 

where the index j is related to time stepping. On the other hand, if one considers 
the Smith method (3.6)-(3.8) for 
S^X + XS = F 

and chooses ^ = j> one obtains the iteration 

(4.12) (I - |s'^)Xj+\l - ^) = (I + *lS'^)Xj(I + llS) + hF . 

In these iterations one may identify the nxn matrix X = [Xj, . . . , xj, Xj € R" and 
the n^ vector 4> = column [Xj^, . . . , xj. If we then identify A^4> with -S^X (i.e., 
Aj = - I®S^) and Aj^* with -XS (i.e. Aj = -S®I), we can immediately see the 
equivalence between (4.11) and (4.12). 



-18- 


5. Summary Remarks and Numerical Examples 

In the prcceeding sections we have presented an algorithm which offers some 
definite advantages in computing directly the feedback gains K for high dimensional 
LQR problems such as those arising in approximating partial or delay differential 
equation control problems. As we shall see with several numerical examples in this 
section, it can substantially outperform standard eigenvector methods on such 
problems. As we have pointed out, a fundamental algebraic operation (in both the 
Chandrasekhar update (2.6), (2.7) and in the reformulated Smith methods (4.10b)) 
involves computation of 

(5.1) L(I - r(A - BK))-i 

where L and K are pxn and mxn matrices respectively. Our algorithm uses the 
Sherman-Morrison-Woodbury formula which can provide significant computational 
savings when m and p are small compared to n. For systems involving sparse 
matrices A (a frequent occurence in many approximation schemes), the needed 
calculations can be carried out quite efficiently. 

We further note that the Chandrasekhar and Newton-Kleinman-Smith 

components as formulated in our algorithm lead to ready estimates between the true 
gain K and the iterates Kj in terms of equation errors in the steps being performed. 

One component (the variable step Smith) of the algorithm is most effectively 
carried out if one possesses some a priori knowledge of bounds on the closed loop 
eigenvalues. If the closed loop eigenvalues lie close to the imaginary axis, then 

convergence in the Smith method can be very slow. Eigen or Schur vector methods 
[27], [30] are less sensitive in this regard. For low order systems, the Schur vector 
approach is more reliable and less expensive computationally than our algorithm. Our 
hybrid algorithm depends critically on a number of choices (e.g., stopping criteria in 
the Newton-Kleinman and Smith components, stepsize sequence (rj) and iteration count 
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sequence {kj} in the variable step component) to be made by the user and the "best" 
choices are heavily problem dependent. Hence one can expect our hybrid algorithm to 
require more experimentation and fine tuning than other more standard methods. 

However, as we shall demonstrate with examples, for the case where n is large 

compared to m and p, it can offer considerable computational savings with no loss in 
accuracy over the methods mentioned above. 

We have tested (and are continuing our efforts in this direction) our hybrid 
algorithm on several numerical examples. We shall report on just two of these here to 
illustrate our findings. All our computations were carried out in double precision on 
an IBM 3081 at Brown University. We gratefully acknowledge the assistance of Yun 

Wang in our carrying out of the extensive computational studies reported for the 

boundary flux control in the diffusion equation problem of Example 2 below. 

Example 1: As one of our examples, we considered an example (Example 6 of [27]) 

which Laub used to test his Schur based methods. The system is the n-dimensional 
system of (1.1) with 


A = 

* o o 

1 0 . 

0 1 

• • 

0 

0 

B = 

0 


0 

• • « 

0 

1 

0 


b 

1 

C = 


0 . . 

. 0 

] 




which leads to an ill conditioned Riccati equation. This problem corresponds to one 
in which n integrators are connected in series with a feedback controller to be 
applied to the nth integrator in order to stabilize the system. Only deviations of Xj 
from the origin are penalized in the cost functional. The true optimal gain is an n 
vector K = (K^ . . . , K") and for this example one can argue that = 1. In [27], 
Laub used his Schur techniques to study this example and reported difficulties with 
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loss of accuracy at a relatively low value of n, n = 21. We carried out runs with 
our hybrid algorithm and obtained quite favorable performance. Some of our findings 
included: 

(a) For n = 40, we used the Chandrasekhar component to integrate to tj = 100 and 
produce an initial estimate KJ = .99041, which when used in the Newton-Kleinman 
(fixed step size r = .5 in the Smith) produces the estimate Kg = 1.0 - in a total of 
2.93 seconds of CPU time. When we used a cruder solution in the Chandrasekhar 
component (tj = 200 but with step size twice that in the first run) to produce Kq = 
.9394, followed by the N - K (r^ = .5, rj = 1.0 in the variable step Smith) we 
obtained Kg = 1.0 -, all in 2.39 seconds. 

(b) For n = 50, we produced kJ = .9224745 at tj = 220 and after the N-K-Smith 
(fixed step r = .5 in the Smith) obtained Kg = 1.0000000003820 in a total of 4.44 
CPU seconds. For the same runs with variable step (r^ = .5, rj = .7) Smith we 
obtained a Kg as above with 3820 replaced by 3817 in a total of 4.31 seconds. 

(c) We compared runs with the Chandrasekhar component only against the Potter 

method for n = 10, 21, 40. Obtaining essentially the same estimates for n = 10 and 
21 (at n = 40, the Potter degenerates numerically to produce useless estimates) we had 
CPU times of = .753 seconds, = .188 seconds, = 1.52 

seconds, POTTj^_ 2 j = 1.22 seconds, CHj^_ 4 q = 4.35 seconds, POTTj^_^q = 6.81 seconds. 

We found for this example that the eigenvector methods are best for small 
n, but as n grows, the Chandrasekhar alone, and, even more so, the hybrid method 
will out perform the eigen-Schur methods in both accuracy and CPU times, A more 
striking demonstration of this behavior will be given in the next example. 

Example 2: We consider the linear quadratic regulator problem: minimize the cost 


functional 
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(5.2) J(u) = j”(|Cz(t)|2 + |u(t)|2)dt 

subject to the partial differential equation 


(5.3) lr^(*’^) = — z(t,x) , X € (0,1) 

0x2 

z(0,x) = <t>(x) 

with boundary conditions 

(5.4) I- z(t,0) = u(t) and |-z{t,l) = 0 

0X 0X 

where c(-) is square integrable on [0,1] and 

Cz(t) = f^c(x) z(t,x) dx 
^0 

We can discretize or approximate (5.3)-(5.4) using the standard Galerkin method [2]; 
i.e. the approximating solution z^(t,x) to (5.3)-(5.4) is given by 

(5.5) zN(t,x) = E W;(t)i .(X) , wj(t) € R\ 

i=0 

where J!. = is the first order spline defined by 


i;N(x) 


N(x - ili— ), li-1. $ X $ L 

N N N 

N(iiL - x), L $ X $ (itll 
N N N 

0 otherwise 


and z^(t,x) satisfies 


(5.6) 


zN(t,x)./.N, 


(x)dx 


- dx - u(t)t/j^0) 

J Q 0X 0X 


for all € Z^=span {i^ ,^ , . . . , } . 

Then, (5.6) leads to the nth order (n = N + 1) ordinary differential equation for 
w^ = co1(Wq, . . . , Wj^) ; 

(5.7) Q^wN(t) = -H^w^(t) - B^u(t) 


where 


22 - 


0 


qn = i 

^ N 


1/3 1/6 0 

1/6 2/3 1/6 

0 


0 

1/6 2/3 1/6 

0 1/6 1/3 


with Qjj = J^Jijdx , 


= N 


1 -1 0 

1 2-1 
0 


0 

-1 2 -1 
0 -1 1 


with ^ ijdx , 

‘J *'odx ‘ dx J 


and 


= col (1 0 ... 0) . 


For computational convenience, we change coordinates (for fixed N) in the 
system (5.7) by x = Q^w^ to obtain the approximate system 
X = -H^(Q^)‘^x - u . 

Thus, in (1.1) we have A = -H^(Q^)‘\ B = -B^ and C = C^(Q^)’^ where is the 
vector with components c|^ = Jj c(x)ij(x)dx, 0 « i ^ N. 

For the problem in this example, the approximating optimal feedback 
operator is given [2] by: 

K^z = k^(x)z(x)dx 

NT N 

where k' (x) = E kjij(x) and K = (kp, . . . , kj^) is the optimal feedback solution in 
i=l 

the problem for (1.1) with A, B, C chosen as indicated above. We note in this case 
that for any N ^ 1, A has only one unstable eigenvalue (zero), (A,B) is stabilizable, 
and (A,C) is detectable. 

For the special case when c(x) = 1, we find C = (1, . . . , 1) e r^^(n+i) 
hence CA = 0. It is thus easy to see that the desired solution (K(t),L(t)) to the 
Chandrasekhar system (1.3) is given by 
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K(t) = k(t)C , L(t) = i(t)C, 
where k,i are scalar functions satisfying 


k = f , k(0) = 0 

i = -jJk , i(0) = 1 


Therefore we find kk + j8jB = 0 so that k^(t) + J!^(t) = 1. We thus find in this case 
that k(t) -♦ 1 as t and hence K = limK(t) = C. For this case, the 

Chandrasekhar system for the infinite dimensional LQR problem (5.2)-(5.4) can also 
be analyzed [17], [36] and exactly the same argument as above shows the optimal 
feedback gain operator is given by 
Kz = l-z(x)dx . 

These analytic solutions can be used to test software packages and approximation 
schemes before more interesting, analytically intractable examples are considered. 
Remark: The form (5.7) of system equations appears frequently in applications. Thus 

the critical computational factor (5.1) can be modified so that one can avoid 
computing A. For example, in this case it has the form 
L(I - r (-HQ-^ - BK))-^ 


(5.8) 


= LQ(Q + rH + rBKQ)-^ 

where Q + rH is a symmetric, tridiagonal, positive matrix. Thus one can readily use 
the Cholesky decomposition algorithm for computing LQ(Q + rH)'^ and combine this 
with the Sherman-Morrison-Woodbury formula (see Remark 3 of Section 2) to 
efficiently compute the critical expression (5.8). 

We carried out extensive computations for this example with c(x) = 1 + x. 
We compared our hybrid method to the Potter algorithm and to use of the 
Chandrasekhar system alone. We have not used the Laub-Schur method on this 
example since we felt comparison with a readily available (to us) Potter package 
would give as a feel for the relative advantages and disadvantages of our scheme 
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compared to eigen and Schur vector based techniques. (Analysis and computational 
experience indicate that the Potter method and the Laub-Schur method are both 
0(N®) with the latter method about twice as fast as the Potter method.) We required, 
whenever feasible, the same level of accuracy in computation of feedback gains and 


compared relative CPU times. 


In studying our hybrid scheme, we tested numerous sets of Smith 
acceleration parameters {rj}, (kj), stopping times t^ in the Chandrasekhar component 
and error stopping criteria in both the Chandrasekhar and Newton-Kleinman-Smith 
components. We summarize some of our findings to date. 

In Table 5.1 we present comparative CPU times for the hybrid scheme vs. 


Potter as we increase N. Recall the corresponding finite dimensional approximation 
scheme has system with dimension n = N+1. In all of the runs reported in Table 5.1, 


the feedback gains for the hybrid and Potter calculations agreed to 9 decimal places 

(CPU Sec.) 


N 

Hvbrid fCPU Sec.I 

Potter 

10 

.17 

.14 

20 

.31 

.81 

30 

.56 

2.45 

40 

.74 

5.49 

50 

.91 

10.71 

60 

1.09 

18.09 

70 

1.26 

27.97 

80 

1.43 

41.56 

100 

1.76 


120 

2.10 


140 

2.45 


160 

2.80 



Table 5.1 

so both scemes provided accurate solutions. In these runs, the hybrid scheme 
calculations used t^ = 2.2 (corresponding to h = .1) with |L(-tf)| “ 10'® in the 
Chandrasekhar component. The Newton-Kleinman component converged after 4 
iterations (i.e. at K^) and we used acceleration steps r^ = 1, rj = 10"^, rg = 10"®, 
r^ = 10‘® . Each Smith iteration was allowed a maximum of kj = 50 per value of Tj 
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although in most cases the iteration satisfied a convergence criterion before this 
maximum was attained. Careful consideration of Table 5.1 reveals that the hybrid 
scheme is clearly 0(N) while the Potter is O(N^); both rates are to be expected from 
our earlier observations about the methods. Note that at N = 80 the hybrid scheme 
is more than 25 times faster than the Potter scheme (with comparable accuracy, of 
course). 

We also ran the hybrid scheme with N = 80 and a number of different 
fixed acceleration values r in the Smith component. The same Chandrasekhar 
component prameters as reported above were used. Table 5.2 contains relative CPU 
times as well as an indication of the N-K iterate for which convergence was 
achieved. 

In Table 5.3 we list some CPU times when different sets of acceleration 
parameters (rj) were used. Again these runs were for N = 80 with the same 


r 



CPU 

: (Sec.) 

Converged 

1 N-K Gain 

5 



5.10 



^6 

1 



6.52 





10'^ 

10-^ 

10-® 

10-* 



6.27 

7.40 

10.10 

12.06 




K4 

^5 

10-® 



10.07 




K4 





Table 

5.2 


r 







CPU (Sec.) 

(10' 

-1 

5 

10 

-2) 




6.25 

(1, 

10- 

-1 

10-2) 




4.88 

(1. 

10' 

-1 

5 

10-2, 

O 

1 

^ CC 

O 

1 

10- 

■®) 

1.98 

(1, 

10' 

-1 

5 

10-^ 

10-^ 10-®, 

10' 

■^) 

1.61 


Table 5.3 

Chandrasekhar solution as above. All of the converged Newton-Kleinman iterates were 
after 6 steps (i.e. Kg). 

Finally, we made runs (for N = 80) to find the best results that the 
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Chandrasekhar algorithm alone (i.e. accurate integration until K(t) K, L(t) -* 0) 
could produce. The best results we were able to achieve yielded an accurate value of 
K for K(-tp with t^ = 3.22 with |L(tj)| “ 10'® obtained in 5.85 CPU seconds. 


Based on our computational findings for the above two examples and our 
experience with several other examples for infinite dimensional systems (e.g., beams 
with tip bodies, etc.), we are quite confident that the hybrid scheme we propose in 
this paper can be profitably used with a number of large scale LQR problems. We 
are currently developing a rather general software package that implements the 
hybrid scheme in a manner so that a broad range of problems can be treated in the 
context of the ideas presented here. 
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